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We investigate the critical behavior of trapped particle systems at the low-temperature superfluid 
transition. In particular, we consider the three-dimensional Bose-Hubbard model in the presence of 
a trapping harmonic potential coupled with the particle density, which is a realistic model of cold 
bosonic atoms in optical lattices. We present a numerical study based on quantum Monte Carlo 
simulations, analyzed in the framework of the trap-size scaling (TSS). 

We show how the critical parameters can be derived from the trap-size dependences of appropriate 
observables, matching them with TSS. This provides a systematic scheme which is supposed to ex- 
actly converge to the critical parameters of the transition in the large trap-size limit. Our numerical 
analysis may provide a guide for experimental investigations of trapped systems at finite-temperature 
and quantum transitions, showing how critical parameters may be determined by looking at the scal- 
ing of the critical modes with respect to the trap size, i.e. by matching the trap-size dependence of 
the experimental data with the expected TSS Ansatz. 



PACS numbers: 05.70.Fh,67.85.-d,67.25.dj,05.30. Jp 
I. INTRODUCTION 

The theory of critical phenomena^ 2 - at phase transi- 
tions generally applies to homogenous systems. How- 
ever, homogeneous conditions are often an ideal limit of 
experimental conditions. Thus, the study of the effects 
of inhomogeneous features is often essential to achieve 
a correct interpretation of the experimental data at the 
transition between different phases, in order to obtain 
reliable estimates of the critical parameters, such as the 
critical temperature, universal critical exponents, etc. 

In this paper we focus on quantum systems of inter- 
acting particles in the presence of an external space- 
dependent potential coupled to the particle density, 
which effectively traps the particles within a limited re- 
gion of space. The presence of a harmonic trap is a com- 
mon feature of the experimental realizations of the Bosc- 
Einstein condensation (BEC) in diluted atomic vapors 3 
and experiments of cold atoms in optical lattices created 
by laser-induced standing waves^, which have provided 
a great opportunity to investigate the interplay between 
quantum and statistical effects in particle systems. 

The critical behavior arising from the formation of 
BEC has been investigated experimentally in a trapped 
atomic system^, observing an increasing correlation 
length compatible with what expected at a continuous 
transition. However, the inhomogeneity due to the trap- 
ping potential drastically changes, even qualitatively, the 
general features of the behavior at a phase transition. For 
example, the correlation functions of the critical modes 
do not develop a diverging length scale in a trap. Never- 
theless, even in the presence of the trap, and in particular 
when the trap gets large, we may still observe a critical 
regime although distorted by the presence of the trap. In 
experiments of trapped particle systems aimed to investi- 
gate their many-body critical behaviors at quantum and 
finite-temperature phase transitions, an accurate deter- 



mination of the critical parameters, such as the critical 
temperature, critical exponents, etc., calls for a quan- 
titative analysis of the trap effects. This issue has been 
much discussed within theoretical and experimental in- 
vestigations, see e.g. Refs. [7rj4li 

Around the transition point, the critical behavior in 
trapped systems is expected to show a power-law scal- 
ing with respect to the trap size, which we call trap-size 
scalin g 17 ! 23 (TSS), controlled by the universality class of 
the phase transition of the homogenous system. TSS has 
some analogies with the standard finite-size scaling (FSS) 
theory for homogenous systems^ 2 - - — , with two main dif- 
ferences: the inhomogeneity due to the space-dependence 
of the external field, and a nontrivial power-law depen- 
dence of the correlation length £ when increasing the trap 
size I at the critical point, i.e. £ ~ I where is the uni- 
versal trap exponent. 

In this paper we show how TSS can be exploited to de- 
termine the critical parameters of trapped particle sys- 
tems, by analyzing the trap-size dependence of observ- 
ables related to the critical modes around the center of 
the trap. The main advantage of this approach is that 
it is supposed to be exact in the large trap-size limit, 
thus providing a systematic scheme to control and im- 
prove the accuracy of the results, without appealing to 
further assumptions and approximations, such as mean- 
field and local-density approximations. As we shall see, 
the method resembles standard FSS techniques, which 
are routinely used to obtain accurate estimates of the 
critical parameters in homogeneous systems, by looking 
at the asymptotic scaling behavior with respect to the 
size of the system, see, e.g., Refs. SEEa 

We investigate this issue at the finite-temperature su- 
perfluid transition of the three-dimensional (3D) Bosc- 
Hubbard (BH) model^, which is particularly relevant 
for cold-atom experiments because it describes bosonic 
atoms in optical lattices^. For this purpose we present 
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FIG. 1: (Color online) Sketch of the T-fJ, phase diagram of 
the 3D BH model ^ in the hard-core U — > oo limit. The 
superfluid transition line starts at (T = 0, fj, = —3), which is 
the location of a quantum transition from the vacuum state to 
the superfluid phase, and ends at another quantum transition 
at (T = 0, fi = 3) (not shown in the figure) between the 
superfluid and Mott phases. The dashed line at a — — 2 shows 
the line along which we numerically investigate the critical 
behavior. 



a numerical analysis based on quantum Monte Carlo 
(QMC) simulations of the 3D BH model with an external 
harmonic potential coupled to the particle density. 

This study may be useful for experimental investiga- 
tions of trapped systems at finite-temperature and quan- 
tum transitions, suggesting some effective recipes to de- 
termine the critical parameters from the scaling of the 
observables related to the critical modes with respect to 
the trap size, i.e. by matching the trap-size dependence 
of the experimental data with the expected TSS Ansatz, 
similarly to experiments probing FSS behavior in 4 He at 
the superfluid transition^. 

The paper is organized as follows. In Sec.|H]we present 
the general features of the TSS in trapped bosonic par- 
ticle systems at the finite-temperature superfluid tran- 
sition driven by the formation of a BEC, such as those 
described by the 3D BH model. Sec. IIIII presents a nu- 
merical study of the superfluid transition of the 3D BH 
model in the hard-core limit, using standard FSS tech- 
niques, which allows us to verify the 3D XY universality 
class of the critical behavior, and to accurately determine 
the critical point. In Sec. IIVI we show how the critical 
parameters can be estimated by a TSS analysis of nu- 
merical QMC data of the trapped BH model. Finally, in 
Sec. [V] we draw our conclusions. App. [A] reports some 
details of the computation of the trap exponent 6 at the 
3D superfluid transition of bosonic systems. 



II. TSS AT THE SUPERFLUID TRANSITION 
OF THE 3D BOSE-HUBBARD MODEL 

A. The phase diagram of the 3D BH model 

Systems of bosonic atoms in optical lattices can be 
modeled by the 3D BH model^^, denned by the Hamil- 



tonian 

#BH = ~$>fo +fc i 6 '-) + (1) 

(ij) 

+ Y^2 n ^ ni ~ ^ ~ vJ2 ni ' 

i i 

where bi is the bosonic operator, m = b\bi is the particle 
density operator, and the sums run over the bonds (ij) 
and the sites i of a cubic lattice. The basic observables 
are related to the particle density 

P(x)^(n x ), (2) 

and the correlation functions of the bosonic field and the 
particle density, 

G b (x,y)^(6+& y ), (3) 
G„(x,y) = (n x n y ) - (n x )(n y ). (4) 

In actual experiments the particle density and its cor- 
relations, such as G n , can be measured by the so-called 
in situ density image techniques, see e.g. Refs. Il4ll49ll50l . 
Some information on the correlation function Gb, and in 
particular the related momentum distribution 

n(k) = £ e ik (*-^G 6 (x, y ), (5) 

can be inferred from the interference patterns of absorp- 
tion images after a time-of-flight period in the large-time 
ballistic regime, see e.g. Ref. |4j. 

The T-fi phase diagram of the 3D BH model presents a 
finite-T transition line separating the normal-fluid phase 
and the low-temperature superfluid phase. The finite-T 
superfluid transition is characterized by the accumula- 
tion of a macroscopic number of bosonic atoms in a single 
quantum state, giving rise to the BEC. The condensate 
wave function naturally provides the complex order pa- 
rameter ip(x) of the phase transition and its relevant U(l) 
symmetry. These global features characterize the 3D XY 
universality class which describes the universal critical 
behavior of a wide class of systems, see, e.g., Ref. |H. 
Numerical studies of the phase diagram of the 3D BH 
model are reported in Refs. I4lll5ll . 

In Fig. Q] we sketch the T-fi phase diagram of the 3D 
BH model in the hard-core U — > oo limit, implying that 
the particle number per site is restricted to the values 
n; = 0,1. The finite-T transition line connects two T = 
quantum critical points at /i — ±3. The T = quantum 
transition at [i — — 3 separates the vacuum state and 
the superfluid phase, while the one at fj, = 3 separates 
the superfluid phase from a Mott phase. In both cases 
the quantum critical behavior is essentially mean field in 
three spatial dimensions^, see also Ref. [52|. 

The critical behavior described by the 3D XY uni- 
versality class is characterized by two relevant param- 
eters t and h, associated with the temperature T, i.e., 
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r ~ T/T c — 1, and the external field h coupled to the 
order parameter. Their renormalization-group (RG) di- 
mensions, y T = \ jv and yh — (5 — fj)/2 respectively, are 
related to the critical exponent v of the correlation length 
and to the exponent r\ describing the power-law decay of 
the two-point function of the order parameter at T c . The 
critical exponents v and r\ are known with great accu- 
racy from theoretical calculations, see, e.g., the results 
reported in Refs. I45ll53l . and experiments at the 4 He su- 
perfluid transition 54 . Recent theoretical estimates of the 
critical exponents ar o 55 ' 56 

v = 0.6717(1), r] = 0.0381(2). (6) 

B. TSS in the presence of the trap 

A common feature of the experiments with cold atoms 4 
is the presence of an external potential V coupled to the 
particle density, which traps the particles within a lim- 
ited space region. In experiments V is usually effectively 
harmonic. We consider a harmonic rotational-invariant 
potential 

V(r) = v 2 r 2 , (7) 

where r = |x| is the distance from the center of the trap, 
which we locate at the origin of the axis, x = 0. This 
trapping force gives rise to a further term in the Hamil- 
tonian: 

#tBH = ff BH + Xy( r *H- ( 8 ) 

i 

Far from the origin the potential V(r) diverges, therefore 
(n.i) vanishes and the particles are trapped. We define the 
trap size by 

I = y/j/v. (9) 

This definition naturally arises 4 ^^ 3 ^ when we consider 
the thermodynamic limit, which is generally defined by 
the limit N, I — > oo keeping N/l 3 fixed (N is the num- 
ber of particles) , and it is equivalent to introducing the 
chemical potential //, as in Eq. (JTJ). In the following, we 
set J = 1 so that I = 1/v. We consider the model at 
fixed chemical potential so we will skip its dependence 
in the following formulas. 

In the presence of a harmonic trap, the large trap- 
size behavior at the transition can be described in the 
framework of TSS. The trapping potential ([7]) coupled 
to the particle density, as in Eq. ((5J, significantly affects 
the critical modes, introducing another length scale ©■ 
Within the TSS framewor k 17 ' 23 , the scaling law of the 
singular part of the free-energy density around the center 
of the trap can be written as 

F sing (x,7» = l-^F{rl- e ,Tl e y^hl eVh )- (10) 

where y T and yh are the RG dimensions reported at the 
end of Sec. Ill Al and 9 is the trap exponent. TSS implies 



that at the critical point (r = 0) the correlation length 
£ of the critical modes is finite, but increases as £ ~ I 6 
with increasing the trap size I. The trap exponent can be 
inferred by a renormalization-group (RG) analysis of the 
perturbation induced by the external trapping potential 
coupled to the particle density. We obtain^ 7 - 

9 = TTto = °- 5732? ( 4 )- ( n ) 

The derivation is outlined in App. [X] 

The TSS equations for the observables and correlation 
functions provide an effective description of the critical 
behavior around the center of the trap, and, in particular, 
of the interplay between the temperature and the confin- 
ing potential. At the superfluid transition and around 
the center of the trap, the one-particle correlation func- 
tion Gb behaves as 

G b (x, y) a r < 1+ "> 9 &(xr fl , yl~° , t1 6 /»), (12) 

and the particle-density correlation G n as 

G„(x, y) a r 2 ^ 9 n (xZ- fl , yr 9 , rl e ' v ), (13) 

where y n is the RG dimension of the density operator 

y n = 3 - l/i/ = 1.5112(2). (14) 

Analogous scaling relations can be inferred for other cor- 
relations. 

In our TSS analyses we consider the trap susceptibility 
Xt defined as 

Xt = ^G b (0,x) (15) 

X 

(we distinguish it from a generic susceptibility because xt 
is the space integral of the correlations with the center of 
the trap), and the trap correlation length £t defined from 
the second moment of Gb(0,x), i.e. 

£ t 2 = -i-£|x| 2 G b (0,x). (16) 
6x* ^ 

According to TSS, they are expected to behave as 

Xt w l^°X{Tl e ' v ), (17) 
6 w l 6 Tl{Tl elv ). (18) 

Note however that any length scale £ extracted from the 
critical modes is expected to show the same TSS as £t. 

The above TSS equations provide the asymptotic de- 
pendence on the trap size I. Scaling corrections are gen- 
erally expected to be 0(Z _we ) where w = 0.785(20) is 
the scaling-correction exponent of the 3D XY universal- 
ity class^S 3 -! 5 ^, thus 

w8 = 0.45(1). (19) 

We remark that the lattice structure of the BH model 
does not play any particular role in the derivation of the 
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TSS formulas, i.e. the microscopic details of the model 
are irrelevant in the TSS limit. Therefore, the above 
TSS equations apply to a wide class of models, i.e. to 
general 3D interacting bosonic systems at the transition 
driven by the BEC, thus also including the atomic system 
experimentally investigated in Ref. [g. 

We also mention that an analogous TSS behavior ap- 
plies to the T = quantum superfluid-to-Mott phase 
transition of the 2D BH model ([1} at fixed integer den- 
sity^, which belongs to the same 3D XY universality 
class 4 ^. 



III. FINITE-SIZE SCALING ANALYSIS 

Before studying the effect of an external space- 
dependent trapping potential, we present a finite-size 
scaling (FSS) analysis of quantum Monte Carlo (QMC) 
simulations of the homogenous BH model (JTJ , with peri- 
odic boundary conditions. In particular, we consider the 
U — » oo hard-core limit of the BH model at fixed chem- 
ical potential fj, = — 2, and vary the temperature along 
the dashed line sketched in Fig. [TJ The QMC simula- 
tions are performed using the stochastic series expansion 
algorithm with the directed operator-loop techniqu e) ' 5S . 



A. Observables and their FSS 

We compute the particle density p, the one-particle 
correlation function Gh(x,y), 59 cf. Eq. ([3]), and the 
particle-density correlation G„(x, y), cf. Eq. Due 
to translation invariance, they only depend on the differ- 
ence of the arguments, i.e. G#(x, y) = G#(x — y). 

We consider the susceptibility 



(20) 



which is the zero-momentum component of the Fourier 
transform of Gb, 



C t (k) = ^e ta G ft (x). 



(21) 



Note that, due to translation invariance, G&(k) coincides 
with the so-called momentum distribution defined by the 
double sum ([3]), apart from a volume factor. When the 
system has periodic boundary conditions, the second- 
moment correlation length £ is conveniently defined by 



density in particle system o 42 ' 62 . In QMC simulations us- 
ing the stochastic series expansion algorithm, T is ob- 
tained from the linear winding number Wi along the i th 
direction, 



T = 



NT - N~ 



(23) 



where N^ and N~ are the number of non-diagonal oper- 
ators which move the particles respectively in the positive 
and negative i th direction. 

FSS predicts the following asymptotic scaling laws of 
the one-particle and particle density correlation functions 
for r = |x| > 0: 

G b (x) «L- 1 -"&(r/L 1 rL 1 /"), (24) 
GnW^L-^Qnir/L^lM"), (25) 



where r 
as 



T/T c — 1. Thus, the susceptibility x behaves 



= L 2 -" 



(26) 



where we have also included the leading 0(L~ U ) scaling 
corrections, and u = 0.785(20) 45 i 53 ^ is the critical ex- 
ponent controlling the leading scaling corrections in the 
3D XY universality class. The dots indicate further scal- 
ing corrections suppressed by higher powers of X/L. The 
scaling functions g and g u are universal apart from a mul- 
tiplicative constant (since \ is 110 1 RG invariant, <?(0) is 
not universal) and a rescaling of the argument. 

We consider the dimensionless RG invariant quantities 



R s = £/L, R r = TL. 



(27) 



According to the FSS theor y 42 i 43 ' 45 , they behave as (see, 
e.g., Ref. [551) 



(28) 



around T c and in the large L limit. / and f u are scal- 
ing functions. In particular the leading one / is uni- 
versal (although it depends on the shape of the volume 
and the choice of the boundary conditions), i.e. it is in- 
dependent of the particular model within the universal- 
ity class, apart from a trivial rescaling of the argument. 
Thus, R* = /(0) is universal. For cubic-shaped lattices 
with periodic boundary conditions, the universal infinite- 
volume limit of i?£ and i?x at T — T c are known with 
great accuracy^ Rt = 0.5924(4) and R% = 0.516(1). 



G 6 (0) - G b (p) 



4sin 2 (p min /2) G 6 (p) 



(22) 



where p = (p m i n ,0, 0), p m i„ = 2ir/L. We also consider 
the so-called helicity modulu o 55 ' 60 T, which is related to 
the spin stiffness in spin models^, and to the superfluid 



B. FSS of QMC data 

We study the FSS of the observables defined above in 
cubic L 3 lattices, up to L = 32 (up to L = 24 for ob- 
servables related to the one-particle correlation function 
Gb), with periodic boundary conditions. 
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FIG. 2: (Color online) QMC data of i? 5 = £/L (bottom) and 
i?r = TL (top) for the 3D homogenous BH model |T} with 
periodic boundary conditions. The vertical dotted line shows 
our final estimate of T c , i.e. T c = 0.7410(1). The horizon- 
tal segments around the crossing point indicate the universal 
asymptotic values^ 7?| = 0.5924(4) and R r = 0.516(1) at T c . 
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FIG. 4: (Color online) i?x (bottom), (middle), and 
X/L 2 ~ v (top) versus tL 1/v with r = T/T c - 1 and T c = 
0.7410, for homogenous 3D BH model with periodic bound- 
ary conditions 
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FIG. 3: (Color online) Data of R r , R s and x/L 2 ~ v at T c 
versus L~" with ui = 0.785. In the case of Ry and R^ we 
also show (by full symbols) their universal L — s> oo limit: 
R T = 0.516(1) and i?| = 0.5924(4). The dotted lines show 
linear fits of the data of Rr and x/L 2 ~ v ■ In the case of R$, 
higher-order scaling corrections appear also important. 



In Fig. [2] we show the QMC results for i?j and Rr ■ 
Their sets of data for different lattice sizes show a clear 
evidence of a crossing point, whose location is expected 
to converge to T c in the large L limit, according to 
Eq. Moreover, the values of and Rr at the cross- 
ing point are consistent with the asymptotic universal 
values i?| and R r reported above. The small deviations 
appear to decrease with increasing the lattice size; they 
are explained by the presence of 0(L~") corrections, see 
Eq. (1251) and also below. An analogous crossing point is 
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shown by the data of the ratio x/L 2 ~ n - 

In order to derive an estimate of T c , we fit the data to 
the Ansatz 

n m 

R = B* + J2 tHT*^" + L~ u bjT j L j / u , (29) 

i=l j=0 

obtained by expanding Eq. (|28|) around r = 0. The best 
estimate of T c is obtained from the data of Rr- Suf- 
ficiently close to T c , for data with \R T /R T — 1| < 0.1 
say, the first terms of the sums [i.e. setting n = 1 and 
to = in Eq. (l29l) ] provide already good fits keeping the 
known universal quantities v, lu and i?^ fixed (in this 
respect the O^L - ^) scaling correction term is necessary 
to achieve fits with acceptable x 2 /d.o.f). For example 
the fit of the data for L > 10 gives T c = 0.74103(1) with 
X 2 /d.o.f « 1.4. We consider 

T c = 0.7410(1) (/x=-2), (30) 

as our final estimate of T c , where the error includes the 
statistical errors of the fits, and takes into account the 
dependence of the results on the choice of the Ansatz and 
the interval of values of T around the transition allowed 
in the fit. Further subleading scaling corrections are con- 
trolled by increasing the minimum value L mm of L of the 
data allowed in the fits. The analysis of the data of R^ 
gives consistent results, but less precise because they are 
apparently affected by larger scaling corrections. 

Fig.Hshows data of % R r and x/L 2 ~ v at T c = 0.741, 
plotted versus L~ u which is the expected order of the 
leading scaling corrections. As expected R^ and Rr con- 
verge to their universal values R T and i?| . The approach 
of Ry and x/L v is approximately linear with respect 
to L~ u , while in the case of R^ also higher-order scal- 
ing corrections appear significant for the available lattice 
sizes. 

Fig. H reports the QMC data of i? T , i? ? and x/L 2 ^ 11 
versus tL 1 ^ with r = T/T c — 1. They show the asymp- 
totic collapse of the data along a universal curve, apart 
from small scaling corrections which get suppressed with 
increasing L. Fig. [5] shows the data of the one-particle 
correlation function Gb at T c , which are consistent with 
the expected asymptotic scaling behavior reported in 
Eq. (El. 

In conclusion, the above FSS analysis of the QMC data 
of the 3D hard-core BH model at /i = — 2 definitely con- 
firms that its superfiuid transition belongs to the 3D 
XY universality class, and provides an accurate deter- 
mination of the (nonuniversal) critical temperature, cf. 
Eq. (HOD. 

C. Finite-size dependence of the particle density 
and its correlators 

The behaviors of the particle density, the compressibil- 
ity, and the particle-density correlators around the tran- 
sition are particularly interesting because they can be 
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FIG. 5: (Color online) L 1+r >G b (n) vs. r/L (where r = |x|) at 
T — T c for homogenous BH systems with periodic boundary 
conditions. The data show the expected scaling behavior (|24|) . 
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FIG. 6: (Color online) The particle density at T c of the 3D 
homogenous BH model. The dashed line shows a linear fit 
of the data to po + cL~ Vn , with po = 0.16187(1) and c = 
0.140(1). 

directly investigated experimentall y 10 ' 14 ' 49 ' 63 ^ 4 , for ex- 
ample by in situ density image techniques. Therefore we 
report a detailed analysis of their data. 

The periodic boundary conditions preserve translation 
invariance, which implies that the particle density p(x) 
is independent of x, thus 

p(x)=p=-l(A0, iV = 5> x . (31) 

X 

Its behavior around the transition point is analogous 
to that of the energy density in spin systems, see e.g. 
Ref. [H, thus 

p^faM+L-y-MrL 1 /") (32) 

at fixed chemical potential, where f a is a nonuniversal 
analytic function of r (and /i), y n is the RG dimension of 
the particle density operator n x , cf. Eq. (1141) . and / s is a 
universal function apart from a factor and a rescaling of 
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for homogenous systems with periodic boundary conditions. 
The dotted line sketches the expected asymptotic behavior 
G„(x) ~ r~ 2yn at small r = |x|. 
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its argument. This behavior is clearly shown by the data 
at T c , which are well approximated by the asymptotic 
formula 



p « po + cL 



(33) 



as shown by Fig. [SJ A linear fit to ([55)1 gives po = 
0.16187(1). 

Fig. [7] reports data of the particle-density correlation 
function at T c . They show the expected scaling behavior, 
obtainable from Eq. (f!?5j) setting r = 0. Note that it 
develops for positive values of G n , while the negative 
data at small distance are pushed toward the origin, not 
contributing to the scaling behavior. The space integral 
of G n gives the compressibility 



(34) 



Its scaling behavior is complicated by the sum around 
x = 0, which gives rise to a nonuniversal analytic contri- 
bution, analogously to the specific heat in 4 He, see e.g. 
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FIG. 9: (Color online) QMC data of the compressibility k, 
at T c . The dashed line shows a linear fit to the predicted 
asymptotic behavior a + bL a '"; in particular, the data for 
L > 8 give a = 0.90(1) and b w -0.80(1) with x 2 / d -°- f - ~ 1-1- 

Ref. HH. Indeed, we expect 

K «. 9a (r)+i Q/ ^ s (rL 1 ^ )j (35) 

where a is the specific heat exponent a = —0.0151(3), 55 
g a is a nonuniversal analytic function of r, and g s is a 
universal function apart from a factor and a rescaling of 
the argument. Notice that, since a < 0, the nonuniversal 
analytic term provides the leading behavior for L — > oo. 
Fig. [8] shows the data of the compressibility. They hint 
at the typical A shape expected in the infinite-volume 
limit, which also characterizes the specific heat at the 
superfluid transition of 4 He, see e.g. Refs. I441I541 . At T c 
they show the asymptotic scaling behavior 



= a + bL a ' v , 



(36) 



see Fig. [9] Linear fits of the available data at T c , up to 
L = 32, gives a « 0.90 and b w -0.80. 



IV. CRITICAL PARAMETERS FROM TSS 

We now consider the 3D BH model in the presence of a 
trapping potential, as described by the Hamiltonian ([5]l. 
We present results of QMC simulations of the 3D hard- 
core BH model at p = — 2 for several values of the trap 
size I, up to I — 14. The trap is centered in the middle 
of a cubic L 3 lattice, with odd L and open boundary 
conditions. More details on the practical implementation 
of QMC simulations of trapped systems can be found in 
Refs. SIM 

The lattice size L is taken sufficiently large to effec- 
tively reproduce the infinite- volume limit, i.e., so that 
the residual finite-size effects can be considered negligi- 
ble compared with the statistical errors. This is checked 
by comparing results at fixed trap size I with increas- 
ing the lattice size L. In particular, after some checks, 
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FIG. 10: (Color online) QMC data of £ t /l e (top) and 
Xt/l {2 ~ v)e (bottom) in the presence of the trap. The verti- 
cal dotted line indicates the critical temperature T c = 0.7410 
obtained from the FSS analysis of Sec. IIII Bl 
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FIG. 11: (Color online) QMC data of £ t /7 e (top) and 
Xt/^ 2 ' 7 ' 9 (bottom) in the presence of the trap vs r/ S//l/ with 
t = T/T c - 1 and T c = 0.741. 



simulations up to I = 10 were generally performed tak- 
ing L = 41 + 1, while for those at / = 14 we considered 
L/l ^ 7. We return to this point below. 



A. TSS analyses of QMC data 

We now show that a TSS analysis of the data for 
trapped systems allows us to determine the critical pa- 
rameters, analogously to the FSS analysis presented in 
the previous section. 

In order to determine T c from TSS, cf. Eqs. (|17|) and 
(fl"8|) . we may exploit the fact that at T = T c (i.e. r = 0), 
the ratios Xt/l 6 ^ 2-71 ^ and £,t/l S become independent of the 
trap size I in the large-/ limit. Therefore, we expect that 
sets of data for different trap sizes cross each other at one 
value of the temperature (apart from scaling corrections), 
providing an estimate of T c , analogously to the FSS anal- 
ysis of the previous section. This is indeed observed in 
Fig. QUI which shows the available data of Xt/l 8 ^ 11 ^ and 
£t/l e versus T. The apparent crossing point of the TSS 
data indicates T c « 0.74. A more accurate estimate is 



achieved by fitting the data to the simple Ansatz 

a + b(T -T c )l 9/u , (37) 

considering data sufficiently close to the crossing point 
to avoid higher powers oItI 8 /". We obtain T c = 0.741(2) 
and T c = 0.742(2) respectively from the data oi^ t /l 6 and 
Xt/^ 2 ~ n ^ 6 ( m both cases from data for I > 8). The error 
takes also into account results from different intervals of 
values of T around T c . The expected O{l~ u0 ) scaling 
corrections, cf. Eq. ([T9|) . appear quite suppressed, at 
least for I > 8. 

Although the available data come from moderately 
large trap sizes, their analysis shows a clear evidence of 
the expected TSS, which allows us to accurately estimate 
the critical temperature T c with a precision of a few per 
mille, in good agreement with the estimate of T c by the 
standard FSS analysis of Sec. ImBl i.e. T c = 0.7410(1). 

In Fig. El we plot the data of of Xt/l 6{2 ~ 9) and £ t /l e 
versus t1 I v . They are consistent with the scaling be- 
havior predicted by Eqs. ([17]) and (|T8")> . approaching a 
universal curve in the large-Z limit. The TSS of the one- 
particle correlation function G(,(0,x) at T c , i.e. 

G fc (0,x) = r( 1+ ")%(A), X = rjl\ (38) 
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FIG. 12: (Color online) The one-particle correlation function 
Gb(0,x) at T c . The dashed line shows the expected small- 
distance behavior G&(0, x) ~ t- _ ' 1+ ' 7 ' ) . 

is nicely reproduced by the data shown in Fig. Q21 At 
small distance r the two-point function Gb is expected to 
show the power-law behavior of the homogenous system, 
i.e. G 6 (0,x) - l/r 1+ ''. 



B. Finite-size effects with the trap 

The finite-size effects in the presence of the trap, i.e. 
when considering the trap within a finite box of size L 
(with open boundary conditions), can be taken into ac- 
count by adding a further dependence on Ll~ e in the 
TSS Ansata^i (JDJ), and in Eqs. (JUJ and (I3J) as well. 
For example the finite-size and trap-size scaling (FTSS) 
of the trap susceptibility xt and the correlation length 
£t, cf. Eqs. (|T5j) and (fTB) . can be written as 

X t ML*-«X(TL x ''',L/l e ), (39) 
it^Ln{rL 1 '\L/l e ). (40) 

The above scaling is confirmed by the data of Xt shown in 
Fig. [T3l obtained by QMC simulations keeping L/l e = 2 
fixed. 

FTSS also implies that at T c the ratio of quantities 
computed in box of size L becomes only function of L/l 
asymptotically, i.e., 

XtihL) 



Xt(l,L-t oo) 



--fx(L/l 9 ), 
k{L/l 6 )- 



(41) 
(42) 



£t(l,L-> oo) 

Their data at T c support this scaling behavior, see 
Fig.IH 

Moreover they tell us that around the transition the 
finite-size effects on xt and £ t get smaller than one per 
mille when L/l e > 7. All data reported in Sec. IIV Al 
which were supposed to correspond to the infinite size 
limit, were obtained by simulations of systems satisfying 
this condition. 



FIG. 13: (Color online) Finite-size scaling of the trap depen- 
dence of the QMC data of x/i 2- " 7 in the presence of the trap, 
from QMC simulations keeping L/l — 2 fixed and using open 
boundary conditions. 



C. Trap-size dependence of the particle density 

Finally, we discuss the trap-size dependence of the par- 
ticle density, which has been considered in the literature 
as a possible probe of critical behavior, due to the exper- 
imental capability of measuring it quite accurately. 

Analogously to the FSS of homogenous systems, the 
scaling behavior of the particle density is more involved, 
because it is dominated by an analytical contribution. In 
the presence of the trap, its behavior is further compli- 
cated by the fact that the particle density depends on 
the distance from the center of the trap. 

To begin with, we consider the particle density at the 
center of the trap, x — 0. Its asymptotic trap-size de- 
pendence is expected to be 



p(0) = 9a(r) + l-^ d g s (Tl e ^) + 



(43) 



where y n 9 = 0.8664(3), g a is a nonuniversal analytical 
function, and g s is a scaling function. Moreover, the 
local-density approximation (LDA), see e.g. Refs. 131371 
suggests that the leading analytical term g a (j) is identi- 
cal to that of Eq. (|32[) for homogenous systems. This is 
supported by the data at T c shown in Fig. [T31 which are 
consistent with the asymptotic formula 



p(0)papo + br 



(44) 



with po equal to the leading constant term of homogenous 
systems, cf. Eq. (3SJ) with p = 0.16187(1). Indeed, a 
linear fit of the data to Eq. (gt]) gives p = 0.1617(3) and 
b = -0.027(1) with x 2 /d.o.f. « 0.4. 

Concerning the space-dependence of the particle den- 
sity, and using rotational invariance, we expect that its 
large trap-size behavior is 

p(x) a f a (r/l, T) + ry» e f s (r/l e ,Tl e /»), (45) 

where f a is again an analytic function. Its analytic de- 
pendence on the ratio r/l is quite natural, because we 
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FIG. 14: (Color online) Finite-size scaling curves of the trap 
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FIG. 16: (Color online) The space dependence of the particle 
density p(x) at T c in the presence of the trap. 
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FIG. 15: (Color online) The particle density at T c and at the 
center of the trap. The full circle along the y-axis shows the 
value po = 0.16187 of the leading asymptotic term obtained 
for homogenous systems, see Fig. [6] The dotted line shows 
a linear fit to po + bl~° Vn , which gives po = 0.1617(3) and 
b= -0.027(1). 



expect it to be a smooth function of 

pi cS (r/l) = f i-V(r) =n-{r/lf, 



(46) 



as suggested by LDA. The asymptotic dependence on r/l 
is shown by the plot of Fig. [16] We expect that the lead- 
ing analytic function f a (r/l,T) is provided by the LDA 
approximation, i.e. by the particle density Ph(fJ-eSiT) of 
the homogenous system in the infinite- volume limit. The 
asymptotic validity of the LDA of the particle density 
was also found at the T = quantum transitions of ID 
and 2D BH models2k2122. 

The above results show that the behavior of the parti- 
cle density around the center of the trap and across the 
transition is quite nontrivial. We finally write it as 

p(x) = p h [n cS {r/l),T] + l-y" 6 f s {r/l B ,t1 6 /»). (47) 

Let us consider the TSS limit at T — T c of this asymp- 
totic behavior, i.e. I — >■ oo keeping X = r/l fixed. Since 
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FIG. 18: (Color online) The density-density correlation 
G n (0,x) at T c . 

in this limit r/l — X/l 1 ^ 6 — > 0, we can expand the ana- 
lytical term, obtaining 

+ r s » 8 / 8 (i) + o(r fl '" +u)9 ), (48) 

where K h (fi) = dp h (p)/dp. Note that the 0(/ -2 ( 1-9 )) 
term cannot be neglected with respect to the scaling 
term, because 

2(1 - 9) < y n 9 < 4(1 - 9). (49) 

Indeed, 2(1-0) = 0.8535(1) &ndy n 6 = 0.8664(3). There- 
fore, in order to determine the universal scaling term, 
we must subtract the terms containing Ph(p) and fc/j(/z). 
They may be evaluated from calculations within the ho- 
mogenous model at p and T fixed, for example consid- 
ering their large-L limit using periodic boundary con- 
ditions. Using the corresponding results at T c for the 
FSS of homogenous systems, see Sec. IIII1 we estimate 
p h (p = -2,T C ) « 0.16187 and n h (p = -2,T C ) « 0.90. 
Then we define the subtracted particle density 

p sub (x) = p(x)-p h (p,T)-l- 2 ^K h (p,T)X 2 

« l-^ e f s (X), X = r/l e . (50) 

Its scaling behavior is nicely confirmed by the corre- 
sponding data plotted in Fig. [TTJ 

Finally in Fig. [18] we show our data for the particle- 
density correlation G n , which vanish at relatively small 
distance, and do not apparently show scaling behaviors, 
likely because the trap size of the available data is still 
too small. Indeed the FSS data of Fig. [7] begin showing 
scaling at relatively large values of the size, essentially 
because the correlation function is significantly nonzero 
only at small distance. 

V. CONCLUSIONS 

We investigate the critical behavior of trapped parti- 
cle systems at the finite-temperature superfluid transi- 



tion driven by the formation of a Bose-Einstein conden- 
sate (BEC). In particular, we consider lattice particle 
systems described by the 3D Bose-Hubbard (BH) model 
©, which is a realistic model of cold bosonic atoms in 
optical lattices^^. We present FSS and TSS analyses 
of numerical QMC simulations of the homogenous and 
trapped 3D BH model in the hard-core U — > oo limit at 
fixed p = —2. 

We show that an accurate study of the critical behav- 
ior, and accurate determinations of the critical parame- 
ters, can be achieved by matching the trap-size depen- 
dence of some appropriate observables with the scaling 
predicted by trap-size scaling (TSS), see e.g. Eqs. (fl"2|) 
and (|18p. The main advantage of this approach is that 
it is supposed to exactly converge to the critical param- 
eters in the large trap-size limit, thus providing a sys- 
tematic scheme to improve the results and control the 
uncertainty, without using further assumptions and ap- 
proximations, such as mean-field and local-density ap- 
proximation (LDA). Although the available QMC data 
are obtained for moderately large trap sizes, our TSS 
analysis allows us to accurately estimate the critical tem- 
perature T c , e.g. T c = 0.741(2) from the analysis of the 
trap correlation length, cf. Eq. (fl~6|). which agrees with 
the critical value T c = 0.7410(1) obtained by a standard 
FSS analysis of QMC simulations of the homogenous 3D 
BH model. 

Our numerical analysis may provide a guide for experi- 
mental investigations of finite-temperature and quantum 
transitions in physical systems that are made inhomoge- 
noues by the presence of an external space-dependent 
field, like the experimentally interesting case of trapped 
atomic systems 3-5 . Indeed, it should show how the criti- 
cal parameters may be determined by looking at the scal- 
ing of the critical modes with respect to the trap size, 
i.e. by matching the trap-size dependence of the exper- 
imental data with the expected TSS Ansatz, similarly 
to experiments probing the finite-size scaling behavior of 
homogenous 4 He systems at the superfluid transition^. 

A few comments are in order concerning the opti- 
mal observables to determine the critical parameters in 
trapped systems. Of course, they are those which are 
mostly determined by the critical modes around the 
center of the trap. In particular, the most convenient 
quantities are those whose leading behavior in the large 
trap-size limit is given by the universal TSS associated 
with the critical modes. In the case of the superfluid 
transition or BEC, the optimal quantities are related to 
the one-particle correlation function of the bosonic field 
around the center x = of the trap, which scales as 
G 6 (0,x) « l-^+^ e g b (rl- e ,Tl e l u ). For example, we con- 
sider the trap susceptibility \ t = Gfc(0, x) and the 
trap correlation length £ t defined from the second mo- 
ment of G & (0,x), which scale as X t ~ l {2 -^ )e X(rl 6 l v ) 
and £t = l e lZ(rl 6 '/"), respectively. Note that any length 
scale £ extracted from the critical modes is expected to 
show the same TSS behavior as £t, and therefore to be 
effective to determine the critical parameters. 
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Since the particle density can be accurately mea- 
sured in experiments of trapped cold atoms, for exam- 
ple by in situ density image technique o 14 ' 49 ' 50 , its be- 
havior may be used as a probe of the critical behavior 
at finite-temperature and quantum transitions, see e.g. 
Refs. [HiSliiilslliil The particle density turns 
out to show a more involved scaling behavior, given by 
Eq. (|45[) . In this case the nontrivial TSS arising from 
the critical modes, which contains the information on 
the critical behavior, does not provide the leading con- 
tribution to the particle density. Indeed, it scales as 
0(l~ Vn6 ) with y n 9 ps 0.866, while the dominant con- 
tribution in the large trap-size limit is given by an an- 
alytical background function f a (r/l,T) related to non- 
critical modes, cf. Eq. (j45|) . Our numerical analysis, 
see Sec. IIV CI shows that such leading analytical back- 
ground is well approximated (actually we conjecture that 
it is exactly given) by the corresponding LDA, i.e. by the 
infinite- volume limit of homogenous systems at chemical 
potential /i c ff = fx — V(r). Thus, TSS provides the lead- 
ing behavior of the deviations from the LDA around the 
transition. Therefore, a careful subtraction is required in 
order to observe the genuine critical term which carries 
information on the critical temperature and the critical 
exponents. As a consequence, in order to infer the crit- 
ical parameters from the particle density, very accurate 
experimental data are required to fit them to Eq. (|45p. 
Alternatively, one may resort to some approximations for 
the analytical background, but this may make the control 
of the real uncertainty more questionable, significantly 
affecting the accuracy of the results. 

The particle-density correlations and compressibility 
are also working optical lattice observables, see e.g. 
Refs. [MfilSi An analogous TSS applies to the 
connected particle-density correlation for example 
G„(0,x) ps l- 29y ^g n {rl- e ,Tl e / v ). Thus information on 
the critical behavior can be achieved by matching the 
experimental data to this TSS Ansatz. However, our nu- 
merical analysis shows that the universal scaling of G n 
at the supcrfiuid transition turns out to appear at rela- 
tively small distance, see Fig. [7l and [TBI which may make 
an accurate determination of its scaling quite hard. 

As already mentioned, a promising quantity is the 
trap susceptibility xt = J2 x Gb(0,x) ~ li 2 -^) 9 with 
(2 - 7])9 fa 1.12, cf. Eqs. (JTSJ and JT7|). However, 
we should note that Xt is n ot proportional to the zero- 
momentum component of the momentum distribution, 
which is given by n(k) = ^ x e lk '^ x_y - ) Gf ) (x, y) even in 
the presence of the trap, and which can be experimen- 
tally related to the interference patterns of absorption 
images after a time-of-flight period in the large-time bal- 
listic regime, see e.g. Ref. |J. Simple considerations show 
that n (k), and in particular its zero-momentum compo- 
nent, is largely dominated by the noncritical regions of 
the trap, while the contribution of the critical modes are 
suppressed, roughly by a total volume factor of the sys- 
tem. Its critical scaling in trapped system in not clear, at 
least in the TSS framework, thus the global momentum 



distribution of the system does not appear promising to 
accurately determine the critical parameters. See, how- 
ever, Refs. [TUlliJiil for a discussion of methods based 
on the measurement of the momentum distribution. 

We remark that the above considerations also apply to 
generic trapped interacting bosonic particle systems at 
the transition driven by BEC, such as the atomic system 
experimentally investigated in Ref. [6j. 

Let us finally mention that a substantial different 
TSS is expected in 2D trapped bosonic systems. The 
finite-temperature superfluid transition of 2D interact- 
ing bosonic systems is described by the Berezinskii- 
Kosterlitz-Thouless (BKT) theory^ 7 -^, which is not as- 
sociated with any spontaneous symmetry breaking and 
emergence of nonvanishing order parameter. Indeed, 2D 
fluids of identical bosons cannot undergo BEC, giving rise 
to a quasi-long range order at sufficiently low tempera- 
ture, characterized by a power-law large-distance decay 
of the one-particle correlation function. Experimental 
evidences of superfluid transitions in trapped quasi-2D 
atomic gases have been reported in Refs. [69h73| . Again, 
the trap gives rise to a substantial distortion of the BKT 
critical behavior of homogenous systems, see e.g. Ref. Hl|. 
The analysis of the BKT renormalizaton-group flow 74 
shows that the asymptotic TSS at the BKT transition 
is characterized by an apparently trivial trap exponent 
9=1, but the asymptotic behaviors show important 
multiplicative logarithms. For example, at T c the trap 
correlation length is expected to increase as £t ~ l t (\nl t ) K 
with k = — 1 in the case of harmonic traps. 
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Appendix A: Derivation of the trap exponent 8 

The trap exponent 9 generally depends on the univer- 
sality class of the transition, on the space dependence of 
the potential, and on the way it couples to the particles. 
In the case of 3D systems of interacting bosonic parti- 
cles at the transition driven by the formation of BEC, 
the value of 9 can be inferred by a renormalization-group 
(RG) analysis of the perturbation Py representing the 
external trapping potential coupled to the particle den- 
sity. Let us consider a generic trapping potential 

V(r) = v p r p (Al) 

where p is a positive even interger number (in the case of 
a harmonic potential p — 2), with a trap size defined as 
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I = J x / p Jv. We follow the field-theoretical approach of 
Refs. [T7II23I that is we consider the 3D <£> 4 quantum field 
theory which represents the 3D XY universality class, see 
e.g. Ref.S 

= J d 3 x [|c^(x)| 2 + r|^(x)| 2 + u|^(x)| 4 ] , (A2) 

where if) is the complex field associated with the order 
parameter, and r, u are coupling constants. Since the 
particle density corresponds to the energy operator |i/>| 2 , 
we can write the perturbation Py as 

P v = J d 3 xV(x)\ip(x)\ 2 . (A3) 

Introducing the RG dimension y v of the constant v of the 
potential (|A1|) . we derive the RG relation 

pyv-P + Vn = 3, (A4) 



where y n = 3 — 1/v is the RG dimension of the den- 
sity/energy operator |?/;| 2 . We eventually obtain 



0=1 = 

Vv 



which gives Eq. (|11[| in the case of a harmonic potential 
with p = 2. Note that Eq. (|A5[) formally gives = 1 in 
the limit p — ¥ 00. This is consistent with the fact that 
the limit p — > 00 corresponds to a homogenous system in 
a spherical box of size I, and open boundary conditions. 
Thus, TSS must reproduce the standard finite-size scal- 
ing of homogenous systems for p — > 00, where the size 
L = I behaves as a scaling field of dimension one in length 
units. 
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